import scipy.io.wavfile as wav
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from scipy import signal, fftpack, interpolate
import numpy as np
import scipy
import wave
import os
import math
from cut_wave import *
import cv2
import librosa.core
from collections import Counter
from vad import *


# nfft = 262144
# nfft = 16384
nfft = 32768
lipOffset = 40
sonicFreq = 20000



def FiltWav(signalWav, N, lowPass, highPass):
    b, a = signal.butter(N=N, Wn=[lowPass, highPass], btype='bandpass')
    filtSignal = signal.filtfilt(b, a, signalWav)  # signalWav为要过滤的信号
    return filtSignal

def MultiVector(A, B):

    C = np.zeros(len(A))
    for i in range(len(A)):
        C[i] = A[i]*B[i]
    return sum(C)

def LMS(xn, M, mu, signalNum):
    en = np.zeros(signalNum)
    W = [0] * M
    for k in range(signalNum)[M - 1:signalNum]:
        x = xn[k - M + 1:k + 1]
        # d = x.mean()
        # print(f'cur {k}')
        # print(W)
        # print(x)
        y = MultiVector(W, x)
        en[k] = xn[k] - y
        W = np.add(W, mu * en[k] * x)

    yn = np.zeros(len(xn))
    for i in range(M):
        yn[i] = xn[i]
    for k in range(len(xn))[M - 1:len(xn)]:
        x = xn[k - M + 1:k + 1]
        yn[k] = MultiVector(W, x)
    return y

def SingleSTFT(filePath, sonicFreq=sonicFreq, lipOffset=lipOffset, nfft=nfft, frameTime=0.15, aheadTime=0.01):
    fs, signalWav = wav.read(filePath)
    print(len(signalWav) / fs)

    lowPass = 2 * (sonicFreq - lipOffset) / fs
    highPass = 2 * (sonicFreq + lipOffset) / fs

    # filtSignal = FiltWav(signalWav=signalWav, N=4, lowPass=lowPass, highPass=highPass)

    signalWav = signalWav / 2**12
    print(f'time max: {np.max(np.abs(signalWav))}')

    filtSignal = signalWav

    # filtSignal = librosa.core.amplitude_to_db(filtSignal)

    #
    #
    # # 归一化
    # print(max(filtSignal))
    filtSignal = filtSignal / max(abs(filtSignal))

    signalNum = len(filtSignal)

    M = 8
    mu = 0.05

    # lmsSignal = LMS(filtSignal, M, mu, signalNum)

    lmsSignal = filtSignal
    nperseg = int(fs * frameTime)
    noverlap = int(fs * (frameTime-aheadTime))
    # f, t, Zxx = signal.stft(lmsSignal, fs, nperseg=nperseg, noverlap=noverlap, nfft=nfft, window=('kaiser', 0.3))
    f, t, Zxx = signal.stft(lmsSignal, fs, nperseg=nperseg, noverlap=noverlap, nfft=nfft)
    # Zxx = Zxx / nfft * 2
    Zxx = CutSTFTBiside(Zxx)
    # Zxx = CalPsd(Zxx)
    Zxx = np.abs(Zxx)

    # print(np.argmax(Zxx, axis=0))

    t = t[0:Zxx.shape[1]]
    realOffset = lipOffset + 10
    down, up = LipGraphPara(sonicFreq, realOffset, fs, nfft)
    f = f[down:up]
    Zxx = Zxx[down:up,:]

    return t, f, Zxx

def SubAdjacent(Zxx, flag=False):
    time = Zxx.shape[1]
    shape = Zxx.shape
    # result = np.zeros((shape[0], shape[1]-1))
    result = np.copy(Zxx[:,:shape[1]-1])
    if flag:
        middle = Zxx.shape[0] // 2
        cutOffset = 3
        for i in range(1, time):
            result[middle-cutOffset:middle+cutOffset, i - 1] = (Zxx[middle-cutOffset:middle+cutOffset, i] - Zxx[middle-cutOffset:middle+cutOffset, i - 1])
        print(result[middle-cutOffset:middle+cutOffset, :])
    else:
        for i in range(1, time):
            result[:, i - 1] = (Zxx[:, i] - Zxx[:, i - 1])
    return result

def LipGraphPara(sonicFreq, lipOffset, fs, nfft, padNum=5):
    span = fs / nfft
    basicPos = int(sonicFreq / fs * nfft)
    freqOffset = int(lipOffset // span) + padNum
    down = basicPos - freqOffset
    up = basicPos + freqOffset
    return down, up

def CalVar(Zxx):
    var = np.var(Zxx, axis=0)
    return var

def GaussWeight(x, sigma):
    w = 1 / ((2 * np.pi) ** 0.5) * np.exp(-x*x/2/sigma/sigma)
    return w

def GaussSmooth(var, lenNum, sigma):
    varNum = len(var)
    varBefore = var[lenNum-1::-1]
    dataAfter = var[sigma-lenNum:]
    varAfter = dataAfter[::-1]
    var = np.concatenate((varBefore, var, varAfter))
    mid = lenNum//2
    coreIndex = [i - mid for i in range(lenNum)]
    coreIndex = np.array(coreIndex)
    w = GaussWeight(coreIndex, sigma)
    result = []
    # for i in range(mid, varNum-mid):
    #     temp = w*var[i-mid, i+]

    for i in range(lenNum, varNum+lenNum):
        temp = np.dot(w,var[i-mid:i+lenNum-mid])
        result.append(temp)
    return np.array(result)

def CalSTFTDiff(filePath, sonicFreq=sonicFreq, lipOffset=lipOffset, fs=44100, nfft=nfft):
    f, t, Zxx = SingleSTFT(filePath)
    down, up = LipGraphPara(sonicFreq, lipOffset, fs, nfft)
    f = f[down:up]
    Zxx = Zxx[down:up, :]
    var = CalVar(Zxx)
    t = t[:Zxx.shape[1]]
    return f, t, var


def Normal(Zxx):
    mean = np.mean(Zxx, axis=1)
    std = np.std(Zxx, axis=1)
    yShape = Zxx.shape[1]
    for i in range(yShape):
        Zxx[:,i] = (Zxx[:,i]-mean)/std

    return Zxx

def CalMean(Zxx):
    mean = np.mean(Zxx, axis=1)
    xShape = Zxx.shape[0]

    for i in range(xShape):
        mask = (Zxx[i, :] > mean[i])
        mask = mask.astype('int64')
        # print(mask.shape)
        # print(Zxx[i,:].shape)
        Zxx[i, :] = Zxx[i, :] * mask
    return Zxx

def Multi2One(Zxx):
    maxNum = np.max(abs(Zxx))
    # print(maxNum)
    scale = 1.0 / maxNum
    Zxx = Zxx * scale
    # print(np.max(Zxx))
    return Zxx

def SumFreq(Zxx):
    time = Zxx.shape[1]
    result = np.zeros((time, ))
    for i in range(time):
        temp = np.sum(Zxx[:, i])
        result[i] = temp
    return result

def GaussSmooth2D(Zxx, kernel, sigma):
    Zxx = cv2.GaussianBlur(Zxx, kernel, sigma)
    return Zxx

def ElimDrift(Zxx):
    time = Zxx.shape[1]
    result = np.zeros(Zxx.shape)
    for i in range(time):
        col = Zxx[:,i]
        minV = np.min(col)
        maxV = np.max(col)
        result[:,i] = (col - minV) / (maxV - minV)
    return result

def GetNoise(filePath):
    data = SingleSTFT(filePath)
    t = data[0]
    f = data[1]
    Zxx = data[2]

    Zxx = np.log10(Zxx)
    time = t.shape[0]
    noise = 0.0
    for i in  range(time):
        mean = np.mean(Zxx[:,i])
        noise += mean
    noise = noise / time
    return noise


def CheckKmeansActive(feature):
    from sklearn.cluster import KMeans
    kmeans = KMeans(n_clusters=2, random_state=0).fit(feature)
    activation = kmeans.labels_
    return activation

def CheckActive(Zxx):
    time = Zxx.shape[1]
    freq = Zxx.shape[0]
    randStart = 1

    noisePercent = 0.85

    activation = np.zeros((time,))

    noise = np.mean(Zxx)
    #################
    # noise = np.mean(Zxx, axis=1)
    # mid = freq // 2
    # offset = 2
    # tempSum = np.sum(noise) - np.sum(noise[mid-1:mid-1+2*offset])
    # noise = tempSum / (freq - 2*offset)
    #################

    if noise < 0:
        noise = noise * (1-noisePercent)
    else:
        noise = noise * (1+noisePercent)

    for i in range(time):
        curFreq = Zxx[:, i]
        #############
        # print(curFreq)
        # tempSum = 0.0
        # for j in range(freq):
        #     if j>=mid-offset+1 and j<=mid+offset:
        #         continue
        #     tempSum += curFreq[j]
        # mean = tempSum / (freq - 2*offset)
        #############
        mean = np.mean(curFreq)
        if mean > noise:
            activation[i] = 1
    return activation

def CheckGradActive(Zxx):
    time = Zxx.shape[1]

    # sumRes = np.var(Zxx, axis=0)
    sumRes = np.sum(Zxx, axis=0)

    lenNum = 10
    sigma = 4
    sumRes = GaussSmooth(sumRes, lenNum, sigma)

    topValue = 2000
    maxValue = np.max(sumRes)
    rate = maxValue / topValue
    sumRes = sumRes / rate

    activation = np.zeros((time,))

    winLen = 8
    maxDiff = 60
    # lowThresh = np.mean(sumRes) * 0.95
    lowThresh = 50

    for i in range(0, time - winLen):
        curWin = sumRes[i:i+winLen]
        curMean = np.mean(curWin)
        curMax = np.max(curWin)
        if curMax - curMean > maxDiff or curMean > lowThresh:
            activation[i:i+winLen] = 1
        # if curMean > lowThresh:
        #     activation[i:i+winLen] = 1

    plt.plot(sumRes)
    plt.plot(activation*1000)
    plt.show()

    return activation

def CheckModelActive(Zxx):
    time = Zxx.shape[1]
    energy = np.sum(Zxx, axis=0)
    activation = np.zeros((time,))

    noiseNum = 10
    varRes = np.var(Zxx[:, 0:noiseNum], axis=0)
    lamdaN = np.mean(varRes)
    print('varRes shape:', varRes.shape)

    divRes = []
    for i in range(time):
        temp = DiveProbe(energy[i], lamdaN)
        divRes.append(temp)
    # topValue = 2000
    # maxValue = np.max(divRes)
    # rate = maxValue / topValue
    # divRes = divRes / rate

    return divRes


def SplitWord(activation, winLen):
    time = activation.shape[0]
    upBound = 2
    deCount = upBound
    start = -1
    index = []
    for i in range(time-winLen):
        window = activation[i:i+winLen]
        judge = (window == 0)
        # print(judge)
        if False not in judge:
            deCount += 1
            if deCount == 2:
                end = min(i + winLen - upBound, time-1)
                index.append((start, end))
                start = -1
            deCount = min(deCount, upBound)
        else:
            if deCount >= upBound:
                start = i
            deCount = 0
    if start != -1:
        index.append((start, time-1))
    return index

#输入的是虚数哦
def CalPsd(Zxx):
    freqLen = Zxx.shape[0]
    conj = Zxx.conjugate()
    res = (Zxx * conj) / freqLen

    return res

def SimpleDopplerShift(Zxx, f):
    avg = np.mean(Zxx)
    sonicPos = f.shape[0] // 2
    diffArr = []
    centerFreq = f[sonicPos]
    maxPos = np.argmax(Zxx, axis=0)
    for i in range(len(maxPos)):
        curFreq = f[maxPos[i]]
        diff = curFreq - centerFreq
        diffArr.append(diff)
    diffArr = np.array(diffArr)
    return diffArr

def OutPic(Zxx, file, saveDir):
    if not os.path.exists(saveDir):
        os.makedirs(saveDir)

    curWord = os.path.splitext(file)[0]
    plt.pcolormesh(Zxx)
    plt.axis('off')
    plt.gca().xaxis.set_major_locator(plt.NullLocator())
    plt.gca().yaxis.set_major_locator(plt.NullLocator())
    plt.subplots_adjust(top=1, bottom=0, right=1, left=0, hspace=0, wspace=0)
    plt.margins(0, 0)
    saveName = f"{saveDir}{os.sep}{curWord}.jpg"
    plt.savefig(saveName)

def SaveAllPic(corpus):
    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]

            Zxx = SubAdjacent(Zxx)
            mean = np.mean(Zxx)
            std = np.std(Zxx)
            Zxx = (Zxx - mean) / std

            # Zxx = GaussSmooth2D(Zxx, (5,5), 2)

            OutPic(Zxx, file, 'save_pic')

def GetEnvelope(signalWav):
    hx = fftpack.hilbert(signalWav)
    hy = np.sqrt(signalWav ** 2 + hx ** 2)
    return hy

def FindPeakOrValley(data):
    dataLen = len(data)
    winLen = 3
    pointIndex = []
    for i in range(data-3):
        curWin = data[i:i+winLen]
        left = curWin[0]
        mid = curWin[1]
        right = curWin[2]
        if (mid > left and mid > right) or (mid < left and mid <right):
            pointIndex.append(i)
    pointValue = data[pointIndex]
    return pointIndex, pointValue

def IsHavePeakOrValley(data, thresh):
    peakIndex, peakValue = FindPeakOrValley(data)
    judge = abs(peakValue) > thresh
    if True in judge:
        return True
    else:
        return False

def GetIndex(envelop, isReverse=False):
    maxValue = max(envelop)
    minValue = min(envelop)
    maxDiff = abs(maxValue - minValue)

    edgeIndex = -1
    if isReverse:
        envelop = envelop[::-1]
    envLen = len(envelop)
    clip = envelop[0:100]
    clipVar = np.var(clip)
    varOffset = 3 * abs(clipVar)
    baseWinLen = 441
    winLen = baseWinLen

    alpha = 0.3

    while (winLen <= 10*baseWinLen and edgeIndex < 0):
        while (alpha >= 0.05 and edgeIndex < 0):
            thresh = alpha * maxDiff
            for i in range(0, envLen - winLen):
                curWin = envelop[i:i + winLen]
                curVar = np.var(curWin)
                curMin = np.min(curWin)
                peakIndex = signal.find_peaks(curWin)[0]
                if len(peakIndex) > 0:
                    peakValue = curWin[peakIndex]
                    judge = (peakValue - curMin) > thresh
                    if True in judge:
                        if curVar > varOffset:
                            edgeIndex = i
                            break
            alpha = alpha - 0.05
        winLen = winLen + baseWinLen

    if edgeIndex < 0:
        edgeIndex = 0

    if isReverse:
        edgeIndex = -(edgeIndex+1)
    return edgeIndex

def CheckEnvActive(envelop):

    lenNum = 10
    sigma = 1
    envelop = GaussSmooth(envelop, lenNum, sigma)

    maxValue = max(envelop)
    minValue = min(envelop)

    start = GetIndex(envelop, False)
    end = GetIndex(envelop, True)

    activation = np.zeros((len(envelop), ))
    activation[start:end] = 1

    plt.plot(envelop)
    plt.plot(activation)
    plt.show()


def ShowGraph(corpus, minTh, winLen, sonicFreq=sonicFreq, lipOffset=lipOffset, fs=44100, nfft=nfft):
    count = 0
    maxCount = 2

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]

            # cutOffset = 8
            # middle = f.shape[0] // 2
            # Zxx = np.concatenate((Zxx[0:middle-cutOffset, :], Zxx[middle+cutOffset:,:]), axis=0)

            # Zxx = np.log10(Zxx)
            # Zxx = GaussSmooth2D(Zxx, 5, 2)

            # index = SplitWord(Zxx, 15, 10)
            # print(index)
            # Zxx = Normal(Zxx)
            # Zxx = CalMean(Zxx)


            Zxx = SubAdjacent(Zxx)
            Zxx = np.abs(Zxx)
            # Zxx = Zxx ** 2
            # Zxx = GaussSmooth2D(Zxx, (5,5), 2)
            # mean = np.mean(Zxx)
            # std = np.std(Zxx)
            # Zxx = (Zxx - mean) / std
            # Zxx[Zxx>1] = 2
            # Zxx[Zxx<=1] = 0

            # Zxx = GaussSmooth2D(Zxx, (5,5), 2)

            # OutPic(Zxx, file, 'save_pic')

            # squareZxx = Zxx ** 2
            # Zxx = squareZxx

            # sinSum = AddFeature(Zxx)

            # Zxx = ElimDrift(Zxx)
            #
            # # 计算方差
            # var = CalVar(Zxx)
            # var = GaussSmooth(var, 20, 4)
            # # index = CondSignal(var, minTh, winLen + 15)
            # #
            # # var = var[index[0]:index[1]]
            # # Zxx = Zxx[:, index[0]:index[1]]
########################
            f = f[0:Zxx.shape[0]]
            t = t[0:Zxx.shape[1]]

            plt.subplot(211+count)
            plt.title(filePath)
            plt.xlabel('Time(s)')
            plt.ylabel('Freqency(Hz)')
            plt.pcolormesh(t, f, Zxx)
            # plt.pcolormesh(t, f, Zxx, cmap='YlOrRd')
            plt.colorbar()
            print(filePath)
            print(f'max: {np.max(Zxx)}, min: {np.min(Zxx)}')
            count += 1
            if count>=maxCount:
                plt.show()
                count = 0

##########################
            # t = t[0:Zxx.shape[1]]
            #
            # plt.subplot(211)
            # plt.title(file)
            # plt.pcolormesh(t, f, Zxx, cmap='YlOrRd')
            # plt.subplot(212)
            # plt.title(file)
            # plt.pcolormesh(t, f, squareZxx, cmap='YlOrRd')
            # plt.show()


def ShowFFT(corpus):
    count = 0
    maxCount = 2

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)

            fs, signalWav = wav.read(filePath)
            lowPass = 2 * (sonicFreq - lipOffset) / fs
            highPass = 2 * (sonicFreq + lipOffset) / fs

            # curbefore = fs // 4
            # curAfter = fs // 4
            # signalWav = signalWav[curAfter:len(signalWav)-curAfter]
            #
            # signalWav = signalWav / max(abs(signalWav))


            t = np.linspace(0, len(signalWav) / fs, len(signalWav), endpoint=False)
            # plt.ylim((-1,1 ))
            # plt.title(f'{file} before')
            # plt.xlabel('time(s)')
            # plt.ylabel('amplitute')
            # plt.plot(t, signalWav)
            # plt.show()
            print(fs)
            # filtSignal = FiltWav(signalWav=signalWav, N=2, lowPass=lowPass, highPass=highPass)
            sigfft = np.fft.fft(signalWav)
            sigfft = abs(sigfft)
            time = np.linspace(0, fs, len(signalWav))

            # plt.ylim((-1, 1))
            # plt.title(f'{file} After')
            plt.xlabel('Time(s)')
            plt.ylabel('Amplitute')
            # 取消边栏
            # plt.axis('off')
            # plt.gca().xaxis.set_major_locator(plt.NullLocator())
            # plt.gca().yaxis.set_major_locator(plt.NullLocator())
            # plt.subplots_adjust(top=1, bottom=0, right=1, left=0, hspace=0, wspace=0)
            # plt.margins(0, 0)
            ########
            plt.plot(time, sigfft)
            plt.show()



def ShowWav(corpus):
    count = 0
    maxCount = 2

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)

            fs, signalWav = wav.read(filePath, False)
            print(len(signalWav))
            lowPass = 2 * (sonicFreq - lipOffset) / fs
            highPass = 2 * (sonicFreq + lipOffset) / fs

            curbefore = fs // 4
            curAfter = fs // 4
            signalWav = signalWav[curAfter:len(signalWav)-curAfter]

            signalWav = signalWav / max(abs(signalWav))


            t = np.linspace(0, len(signalWav) / fs, len(signalWav), endpoint=False)
            # plt.ylim((-1,1 ))
            # plt.title(f'{file} before')
            # plt.xlabel('time(s)')
            # plt.ylabel('amplitute')
            # plt.plot(t, signalWav)
            # plt.show()

            filtSignal = FiltWav(signalWav=signalWav, N=2, lowPass=lowPass, highPass=highPass)
            filtSignal = signalWav
            # envelope = GetEnvelope(filtSignal)

            # plt.ylim((-1, 1))
            # plt.title(f'{file} After')
            plt.tick_params(labelsize=12)
            plt.xlabel('Time(s)', fontsize=16)
            plt.ylabel('Amplitute', fontsize=16)
            # 取消边栏
            # plt.axis('off')
            # plt.gca().xaxis.set_major_locator(plt.NullLocator())
            # plt.gca().yaxis.set_major_locator(plt.NullLocator())
            # plt.subplots_adjust(top=1, bottom=0, right=1, left=0, hspace=0, wspace=0)
            # plt.margins(0, 0)
            ########
            # plt.title(filePath)
            plt.plot(t, filtSignal)
            plt.show()




            # filtSignal = np.abs(filtSignal)
            # # plt.ylim((-1, 1))
            # plt.title(f'{file} After')
            # plt.xlabel('time(s)')
            # plt.ylabel('amplitute')
            # plt.plot(t, filtSignal)
            # plt.show()

            # cutlittle = fs // 10
            # t = t[cutlittle:len(envelope)-cutlittle]
            # envelope = envelope[cutlittle:len(envelope)-cutlittle]
            # clipIndex = np.arange(0, len(envelope), 10)
            # envelope = envelope[clipIndex]
            # activation = CheckEnvActive(envelope)
            # print(len(envelope))
            # plt.plot(t, envelope)
            # plt.show()

def ShowFreq(corpus):
    count = 0
    maxCount = 2

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]

            # plt.title(filePath)
            # plt.xlabel('time(s)')
            # plt.ylabel('freqency(hz)')
            # plt.pcolormesh(t, f, Zxx)
            # plt.show()

            print(f'freq max: {np.max(np.abs(Zxx))}')

            Zxx = SubAdjacent(Zxx)
            Zxx = np.abs(Zxx)
            t = t[0:Zxx.shape[1]]

            # Zxx = np.log10(Zxx)

            print(Zxx.shape)
            # plt.title(filePath)
            plt.tick_params(labelsize=12)
            plt.xlabel('Time(s)', fontsize=16)
            plt.ylabel('Freqency(Hz)', fontsize=16)
            # 取消边栏
            # plt.axis('off')
            # plt.gca().xaxis.set_major_locator(plt.NullLocator())
            # plt.gca().yaxis.set_major_locator(plt.NullLocator())
            # plt.subplots_adjust(top=1, bottom=0, right=1, left=0, hspace=0, wspace=0)
            # plt.margins(0, 0)
            ########
            # plt.title(filePath)
            plt.pcolormesh(t, f, Zxx)
            plt.show()

# def ShowFFT(corpus):
#
#     for root, dirs, files in os.walk(corpus):
#         for file in files:
#             filePath = os.path.join(corpus, file)
#
#             fs, signalWav = wav.read(filePath)
#             lowPass = 2 * (sonicFreq - lipOffset) / fs
#             highPass = 2 * (sonicFreq + lipOffset) / fs
#
#             curbefore = fs // 4
#             curAfter = fs // 4
#             signalWav = signalWav[curAfter:len(signalWav) - curAfter]
#
#             signalWav = signalWav / max(abs(signalWav))
#
#             t = np.linspace(0, len(signalWav) / fs, len(signalWav), endpoint=False)
#
#             filtSignal = FiltWav(signalWav=signalWav, N=4, lowPass=lowPass, highPass=highPass)
#
#             fftResult = np.fft.fft(filtSignal)
#             fftResult = np.abs(fftResult)
#             fftResult = fftResult[0:len(fftResult)//2]
#             index = np.argmax(fftResult)
#             print(index)
#             fftResult = np.log10(fftResult)
#             # fftResult[index-2:index-2] = 0
#
#
#             plt.title(f'{file} fft')
#             plt.xlabel('frequency(Hz)')
#             plt.ylabel('amplitute')
#             plt.plot(fftResult)
#             plt.show()

def ShowDiffSplit(corpus):
    count = 0
    maxCount = 2

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]
            plt.title(file)
            Zxx = SubAdjacent(Zxx)
            squareTemp = Zxx ** 2
            sumRes = CheckGradActive(squareTemp)
            divRes = CheckModelActive(squareTemp)
            plt.plot(sumRes)
            plt.plot(divRes)
            plt.show()

def ShowWavFreq(corpus, minTh, winLen, sonicFreq=sonicFreq, lipOffset=lipOffset, fs=44100, nfft=nfft):

    count = 0
    maxCount = 2

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            fs, signalWav = wav.read(filePath)
            timeAxe = np.linspace(0, len(signalWav)/fs, len(signalWav), endpoint=False)
            # plt.subplot(211)
            # plt.title(file)
            # plt.xlabel('time(s)')
            # plt.ylabel('amplitude')
            # plt.plot(timeAxe, signalWav)

            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]
            Zxx = CalPsd(Zxx)
            # Zxx = np.log10(Zxx)
            # Zxx = SubAdjacent(Zxx)

            # cutOffset = 3
            # middle = f.shape[0] // 2
            # Zxx = np.concatenate((Zxx[0:middle-cutOffset, :], Zxx[middle+cutOffset:,:]), axis=0)

            f = f[0:Zxx.shape[0]]
            t = t[0:Zxx.shape[1]]

            # plt.subplot(212)
            plt.title(file)
            plt.xlabel('time(s)')
            plt.ylabel('freqency(hz)')
            plt.pcolormesh(t, f, Zxx, cmap='YlOrRd')
            plt.show()

def ShowMicroMark(corpus, winLen, startFreq, endFreq, interval, lipOffset=lipOffset, nfft=nfft):
    lineWidth = 4
    # syLen = 8

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)

            fs, signalWav = wav.read(filePath)
            startCut = int(0.4 * fs)
            endCut = int(0.12 * fs)
            signalWav = signalWav[startCut:len(signalWav) - endCut]

            for sonicFreq in range(17650, 18000+1, interval):
                data = SingleSTFT(filePath,
                                  sonicFreq=sonicFreq,
                                  lipOffset=lipOffset,
                                  nfft=nfft)
                t = data[0]
                f = data[1]
                Zxx = data[2]

                plt.title(file + " "+ str(sonicFreq))
                Zxx = SubAdjacent(Zxx)
                squareTemp = Zxx ** 2
                activation = CheckGradActive(squareTemp)
                print(f'activation len: {len(activation)}')
                # divRes = CheckModelActive(squareTemp)
                # plt.plot(divRes)
                # plt.show()

                index = SplitWord(activation, winLen)

                # for case in index:
                #     clip = activation[case[0]:case[1]]
                #     cal = Counter(clip)
                #     if cal[1] < syLen:
                #         activation[case[0]:case[1]] = 0
                #         index.remove(case)

                print(filePath)
                print(index)

                t = t[0:Zxx.shape[1]]
                f = f[0:Zxx.shape[0]]
                # oriZxx = oriZxx[0:f.shape[0], 0:t.shape[0]]

                # timeTime = np.linspace(0, len(signalWav) / fs, len(signalWav), endpoint=False)
                # plt.subplot(211)
                # plt.title(file)
                # plt.xlabel('time(s)')
                # plt.ylabel('amplitude')
                # plt.plot(timeTime, signalWav)
                #
                # plt.subplot(212)
                # plt.title(file)
                # plt.xlabel('time(s)')
                # plt.ylabel('freqency(hz)')
                # plt.pcolormesh(t, f, Zxx, cmap='YlOrRd')
                # line = np.zeros((Zxx.shape[0], 2))
                # for i in range(t.shape[0] - 1):
                #     if activation[i] == 1:
                #         plt.pcolormesh(t[i:i + 2], f, line)
                # plt.show()
                #########

def ShowMark(corpus, minTh, winLen, sonicFreq=sonicFreq, lipOffset=lipOffset, fs=44100, nfft=nfft):

    lineWidth = 4
    syLen = 8

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)

            fs, signalWav = wav.read(filePath)
            startCut = int(0.4 * fs)
            endCut = int(0.12 * fs)
            signalWav = signalWav[startCut:len(signalWav)-endCut]
            lowPass = 2 * (sonicFreq - lipOffset) / fs
            highPass = 2 * (sonicFreq + lipOffset) / fs
            signalWav = FiltWav(signalWav=signalWav, N=4, lowPass=lowPass, highPass=highPass)

            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]

            plt.title(file)
            # oriZxx = np.log10(Zxx)
            Zxx = SubAdjacent(Zxx)
            squareTemp = Zxx ** 2
            activation = CheckGradActive(squareTemp)

            index = SplitWord(activation, winLen)

            for case in index:
                clip = activation[case[0]:case[1]]
                cal = Counter(clip)
                if cal[1] < syLen:
                    activation[case[0]:case[1]] = 0
                    index.remove(case)

            print(filePath)
            print(index)

            # print(filePath)
            # print(index)
            # print(activation)
            # print(activation.shape)
            ########
            t = t[0:Zxx.shape[1]]
            f = f[0:Zxx.shape[0]]
            # oriZxx = oriZxx[0:f.shape[0], 0:t.shape[0]]

            # plt.subplot(211)
            # plt.title(file)
            # plt.xlabel('time(s)')
            # plt.ylabel('freqency(hz)')
            # plt.pcolormesh(t, f, Zxx, cmap='YlOrRd')

            timeTime = np.linspace(0, len(signalWav)/fs, len(signalWav), endpoint=False)
            plt.subplot(211)
            plt.title(file)
            plt.xlabel('time(s)')
            plt.ylabel('amplitude')
            plt.plot(timeTime, signalWav)

            plt.subplot(212)
            plt.title(file)
            plt.xlabel('time(s)')
            plt.ylabel('freqency(hz)')
            plt.pcolormesh(t, f, Zxx, cmap='YlOrRd')
            line = np.zeros((Zxx.shape[0], 2))
            for i in range(t.shape[0]-1):
                if activation[i] == 1:
                    plt.pcolormesh(t[i:i+2], f, line)
            plt.show()
            #########


def ShowMiddle(corpus, minTh, winLen, sonicFreq=sonicFreq, lipOffset=lipOffset, fs=44100, nfft=nfft):

    count = 0
    maxCount = 2

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]

            Zxx = SubAdjacent(Zxx)

            f = f[0:Zxx.shape[0]]
            t = t[0:Zxx.shape[1]]

            middle = f.shape[0] // 2

            plt.subplot(211+count)
            plt.title(file)
            plt.xlabel('time(s)')
            plt.ylabel('freqency(hz)')
            plt.plot(t, Zxx[middle, :])
            count += 1
            if count>=maxCount:
                plt.show()
                count = 0


def ShowVar(corpus, minTh, winLen, sonicFreq=sonicFreq, lipOffset=lipOffset, fs=44100, nfft=nfft, lenNum=20, sigma=4):

    count = 0

    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]

            Zxx = SubAdjacent(Zxx)
            t = t[0:Zxx.shape[1]]
            # Zxx = GaussSmooth2D(Zxx, (5,5), 2)

            var = CalVar(Zxx)
            # var = SumFreq(Zxx)
            var = GaussSmooth(var, lenNum, sigma)

            filtIndex = CondSignal(var, minTh=minTh, winLen=winLen)

            plt.subplot(211+count)
            plt.title(file)
            plt.plot(var)
            plt.subplot(212)
            plt.title(file)
            plt.plot(var[filtIndex[0]:filtIndex[1]])
            plt.show()

def ShowSmooth2D(corpus, minTh, winLen, sonicFreq=sonicFreq, lipOffset=lipOffset, fs=44100, nfft=nfft, lenNum=5, sigma=2):

    count = 0
    lenNum = (lenNum, lenNum)
    for root, dirs, files in os.walk(corpus):
        for file in files:
            filePath = os.path.join(corpus, file)
            data = SingleSTFT(filePath)
            t = data[0]
            f = data[1]
            Zxx = data[2]

            Zxx = SubAdjacent(Zxx)
            # Zxx = np.log10(Zxx)
            print(Zxx.dtype)
            t = t[0:Zxx.shape[1]]


            Zxx = GaussSmooth2D(Zxx, (5 ,5), 2)
            # minVal = np.min(Zxx)
            # Zxx = Zxx + np.abs(minVal)
            # avgVal = np.mean(Zxx[0, 0:5])
            # mask = abs(Zxx - avgVal) > 0.5
            # print(mask.shape)
            # mask = mask.astype(np.int)
            # print(mask)
            # Zxx = Zxx * mask
            # mean = np.mean(Zxx)
            # std = np.std(Zxx)
            # Zxx = (Zxx - mean) / (std + 0.0001)
            # minVal = np.min(Zxx)
            # maxVal = np.max(Zxx)
            # Zxx = (Zxx - minVal) / (maxVal - minVal)

            # print(f'{file} maxValue: {max(Zxx[:,5])}')
            plt.title(file)
            plt.pcolormesh(t, f, Zxx, cmap='YlOrRd')
            plt.colorbar()
            plt.show()

def getFFTFeature(filtSignal, fs, sonicFreq, lipOffset):
    sigLen = len(filtSignal)
    maxLen = fs * 2
    if (sigLen > maxLen):
        print('current sigLen surpass maxLen, which is', sigLen)
        assert sigLen > maxLen
    fftResult = np.fft.fft(filtSignal)
    # 数据归一化
    fftResult = fftResult / (sigLen / 2)
    # 插值，将所有数据都对齐
    oldx = np.linspace(0, fs, sigLen)
    newx = np.linspace(0, fs, maxLen)
    f = interpolate.interp1d(oldx, fftResult, kind='quadratic')
    newfft = f(newx)
    # 对fft进行log处理，拉高低值
    newfft = np.log10(newfft)
    # 找到超声波高峰位置，取±40Hz的区域
    peakPos = int((sonicFreq / fs) * maxLen)
    perBin = fs / maxLen
    offPoint = math.ceil(lipOffset / perBin)
    lipfft = newfft[peakPos-offPoint:peakPos+offPoint+1]
    return lipfft


def getStaticFeatureFromUltra(filePath, sonicFreq, lipOffset):
    fs, signalWav = wav.read(filePath)

    lowPass = 2 * (sonicFreq - lipOffset) / fs
    highPass = 2 * (sonicFreq + lipOffset) / fs

    filtSignal = FiltWav(signalWav=signalWav, N=4, lowPass=lowPass, highPass=highPass)

    filtSignal = filtSignal / (2^16)

    envelope = GetEnvelope(filtSignal)
    maxValue = np.max(envelope)
    minValue = np.min(envelope)
    mean = np.mean(envelope)
    std = np.std(envelope)

    fftFeature = getFFTFeature(filtSignal, fs, sonicFreq, lipOffset)
    print(fftFeature.shape)
